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Abstract 

I present a highly efficient method for evolving parton distributions in perturbative QCD. The 
method allows evolving the parton distribution functions according to any of the commonly-used 
truncations of the evolution equations (which differ in their treatment of higher-order terms). I 
also give formulae for computing crossing functions within the method. 
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1. Introduction 


The parton distribution functions of the nucleon are fundamental ingredients in the descrip¬ 
tion of short-distance scattering experiments involving hadronic initial states, whether these be in 
hadron-hadron collisions, or in deeply inelastic scattering. Indeed, along with the running coupling 
CKs; they are the only ingredients needed from outside perturbation theory for a complete description 
of such processes ignoring subleading power corrections. 

Both the parton distributions and the running coupling depend on the momentum scale at 
which they are evaluated. As is well known, the variation in each of these quantities as one moves 
from one momentum scale to another is described by an evolution equation that can be computed 
perturbatively. It is only the values of the distributions and of the coupling at a fixed scale which 
are required inputs from outside perturbation theory. 

Lattice gauge theorists may calculate these quantities someday, but for the moment they must 
be extracted from experiments. The modern approaches [1,2] involve global fits to all available 
experiments. The experiments in fact involve different scale arguments to the distribution functions, 
but as these are related by a perturbative evolution equation, we can regard the fits as determining 
the distributions at a certain fixed scale Qq. 

Present theoretical abilities allow extensive calculations at next-to-leading order (NLO), and 
there are growing indications that next-to-next-to-leading order (NNLO) calculations are required 
if we are to determine ag to 1%. To the required precision, the evolution equation for the running 
coupling can be solved in closed form. This is not true, however, of the evolution equation for the 
parton distributions; these must be evolved numerically from the original fixed scale. 

Workers have used two main approaches to evolve the distribution functions numerically: 
direct integration of the evolution equations, and an approach based on Mellin transformations. 
This paper presents an improvement to the techniques heretofore used, within the framework of 
the latter approach. 


2. Evolution Equations 


A parton distribution function f{x,Q‘^) obeys an evolution equation [3] in momentum. 




dQ^ 


= P{x,Q ) (g) f{x,Q ) 


( 2 . 1 ) 


where P{x,Q‘^) is the Altarelli-Parisi kernel, and (g) denotes a convolution, 

dy dz 6{x - yz)A{y)B{z), 
Jo 


A{x) (g) B{x) = 


( 2 . 2 ) 
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(For a different early discussion, see ref. [4].) For the quark non-singlet distributions, each of the 
quantities in eqn. (2.1) is a scalar; in the singlet sector, / is a two-component vector containing the 
quark singlet distribution and the gluon distribution, and P(x, Q'^) is a 2 x 2 matrix. 

The kernel P{x,Q^) has a perturbative expansion, 

P{x,Q^) = as{Q^)Po{x) + al{Q^)Pi{x) + 0{a^,), (2.3) 


in which as{Q‘^) = Q;s(Q^)/(47r) is a rescaled version of the usual running coupling. As mentioned 
in the introduction, one approach to evolving the parton distribution functions /, known as the 
rc-space method, is to integrate the perturbative approximation to eqn. (2.1) directly. I will return 
to the issues surrounding the choice of method in section 8. For the moment, let us proceed by 
Mellin transforming the evolution equation; with Mellin moments of a function h{x) defined via 


= f dx x^ ^h{x ), 

Jo 


the evolution equation, truncated to next-to-leading order, becomes 

■)2 


= [asmP^ + aUQ^)Pj] riQ^ 


(2.4) 


(2.5) 


since the Mellin transformation turns convolutions into products. (That is reason for using it.) It 
is convenient to change variables, and to use Ug as the evolution variable. We can do this using the 
beta function, 


fi{as) = 




obtaining 


= -Poa'i - Pial + 0{at) , 

r- 


Sin 

dr _ [Po" + a.Pf] 


( 2 . 6 ) 


(2.7) 


dag ag[PQ+ PiQg^ 

The usual approach [5] expands the right-hand side, in keeping with a perturbative treatment of 
the equation, upon which we obtain 


df" 

dag 


1 


PoQg 




r 


( 2 . 8 ) 


The boundary condition is /^(Qq) ~ fo- Conventionally one proceeds by introducing an 
evolution operator E such that 


nQ^) = E^{ag{Q^),ag{Ql))r{Ql). 


(2.9) 


In the non-singlet case, E is just a scalar function of z; in the singlet case, it is a 2 x 2 matrix. The 
evolution operator satisfies the same equation (2.8) we have been considering above, 

dE^ 1 


dag 


Poo-g 


Pn -|- ag 


pz Pi pz 

■ Wo^^ 


E^ 


( 2 . 10 ) 
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but has a boundary condition, E^{ao,ao) = 1, that renders it an entirely perturbative object. 

Solving the evolution equation is straightforward in the non-singlet case; expanding in Og 
beyond the leading terms, we obtain [6,7] 


Erj (Ug, Uq) 



Os — flO . 



( 2 . 11 ) 


where r/ = ±1 correspond respectively to the combinations g — g and {(u + u) — {d + d ),...}, and 
where Pl{r]) = Piiv) ~ singlet case, we need hrst to dehne two matrices which project 

onto the eigenvectors of the leading-order AP coefficient matrix. 


= ± 


(V - Al) 


{p; - a9 


in which 

A"± = 

are the eigenvalues of Pq . 
sector as [6,7], 


pz p2 + 

^0,gg + ^0,qq ^ 


'J{P«.„ - Pl,P + iP, 


pz 

0,qg-^0,gq 


( 2 . 12 ) 


With these definitions, we can write the evolution operator in the singlet 


E (cts, ^o) 


jyz _ ^0 pz pz 


Po 


f3o — + Xt y 

+ (+--)• 


V«o/ 


(A1-A;)//3o 


-ao PlP^Pi 


-Al//3o 


(2.13) 


With the normalizations used here and Uf the number of active flavors, the beta function [8] 
coefficients are 

2 38 

/3o = 11 - -Uf, (3i = 102 - —Uf ; (2.14) 

the Mellin moments of the Altarelli-Parisi function may be found in ref. [9], noting that Pf = 
The analytic continuations of these moments are given by Gliick, Reya, and Vogt [7]; I 
discuss the analytic continuation of a certain term involving the dilogarithm in an appendix. 

Of course, for use in perturbative cross section calculations, we want the evolved distributions 
as functions of the parton momentum fraction x, not their Mellin transforms as functions of the 
conjugate variable z. To obtain the evolved parton distribution in x-space, we must invert the 
Mellin transform, 

f{x,Q^) = ^J^dzx-^nQ^), (2.15) 

where the contour C runs parallel to the imaginary axis, to the right of the rightmost pole in the 
integrand. 
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Up to this point, everything I have written is completely standard. I now wish to investigate 
how one may deform the contour C in order to obtain an efficient method for evaluating the integral. 
Note that the integrand has poles only along the real axis; denote the rightmost pole by c^. (It is 
typically around 0.5 for the non-singlet case, and 1.3 for the singlet case.) 

One can in principle perform the contour integral along the textbook contour z = c + iy 
(—00 < y < 00 and c > c^), as was done in the original work of Gliick and Reya [10]. One then 
obtains the expression 


X ^ 


/ dy {x-^yr+^y{Q‘^) + x^yr-^y{Q‘^)) 

Jo 

J.-C fOO 

= — / dy Re [x-^yr^^y{Q^) + x^yr-^y{Q^)] 
^ Jo 


(2.16) 


The large-y behavior of the integrand is determined by the behavior of the initial distribution 
as x ^ 1. Initial distributions of the form x“(l — x)^ give rise to a power decay, / ~ y~^~^ as 
y ^ 00. Because the contour is parallel to the imaginary axis, the x~^ factor is a purely oscillating; 
it does not damp the integrand as y ^ 00 . On the other hand, because the integrand has no poles 
off the real axis, and because the integral along a half-circle at 00 in the left half-plane (or along part 
of this half-circle) vanishes, we can freely deform the contour into the left-half plane, so long as we 
stay away from the real axis. Were we to choose a contour such that z has an increasingly negative 
real part as it heads off to infinity, the x~^ factor would contribute an exponential suppression, 
improving the convergence of the integral. (Recall that the base point x lies between 0 and 1.) 

This motivates the choice of contour in ref. [7], where it is taken to contain two line segments 
running diagonally into the left-half plane from a point c on the real axis to infinity. 

However, we might ask: why this contour? Or, phrased differently: how should we choose a 
contour? I address this question in the next section. 


3. Choosing a Contour 

The most obvious answer is that we should choose the contour of steepest descent; this choice 
will yield an integral that converges most efficiently and (one therefore hopes) can be evaluated 
numerically with fewest function evaluations for a fixed desired accuracy. One might fear that 
finding the contour of steepest descent requires a complicated computation. For our problem, 
however, it will turn out that there is a simple but very good approximation to the desired contour. 

Now, we don’t want to find a contour for each value of Q^, especially if finding such a contour 
involves a significant number of function evaluations; this would defeat the whole purpose of the 
exercise! For use in a program evaluating a cross section, it is convenient to set up a grid of x and 
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points {a la MRS [1]), and interpolate between them. Instead of finding a contour for each 
pair, use the following strategy: for each grid x value, find the contour of steepest descent 
for Qq, and then use it for evolution to all values. Since the evolution is relatively slow (it is only 
in In In 5 after all), the contour at Qq should be fairly close to the contour of steepest descent 
at Q^. From a programming point of view, this approach allows all the contours, points along 
them, and the associated anomalous dimensions to be computed as setup code; only the evolution 
operators at the given points will need to be evaluated anew for each to which we wish to evolve 
the distributions. 

Let us then first consider the inverse Mellin transform of /^(Qq)- (Yes, we already know the 
answer — it is just f{x, Qq) — but that is not the point.) There are certain technical complications 
in the singlet sector, which I will discuss section 5, so let us restrict attention here to the non-singlet 
distributions. 

The first step in determining the contour is to find the minimum of the integrand along the 
real axis. As the starting distribution is a positive function, its Mellin transform will be positive 
along the real axis to the right of its rightmost pole c^. Furthermore, the starting distribution is 
not infinitely concentrated at x = 1, so its Mellin transform will decrease as z ^ oo. On the other 
hand, since x < 1, x~^ increases exponentially as z ^ oo. The product of the two must therefore 
have a minimum somewhere, and the value of the product there will be positive. (If there happen 
to be multiple minima, pick the one closest to the rightmost pole.) Call this point cq. 

The integrand is analytic as a function of z. The minimum of the function along the real axis 
is therefore actually a saddle point of the function in the complex plane, and thus the place to start 
tracing out our desired contour. Furthermore (again because of analyticity), the desired contour is 
also a contour of stationary phase, so that the integrand will be real along it. (I will also assume 
that the rightmost singularity of the Mellin transform of the initial distribution is to the right of 
the rightmost pole in the anomalous dimensions; this assumption, which holds for realistic parton 
distributions, ensures that the minimum is only slightly different for /^(Q^) than for f^iQo), and 
is in any event necessary if we are to use a contour as determined below for integrating /^(Q^).) 

Define F{z) = x~^(Q q)', our inverse Mellin transform then has the form 


I = 


2TTi 


’Cs + i-Cs) 


dzx ^f\Ql) = ^ f [dz F{z) - dz F{z)] 

^ C s 


(3.1) 


TT 


Re [—1 dz F{z)] 




where Cs is the part of the contour of steepest descent running upwards from cq. 

Let us examine the contour near the point cq. A generic saddle point will have a non-vanishing 
second derivative (and all the relevant minima for conventional choices of /(x, Qq) are indeed 
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generic), so the contour will start out with a tangent parallel to the imaginary axis: that is the 
direction of steepest descent of 

Fiz) ~ F(co) + (z - co)2 + • • • , (3.2) 

which is perhaps easiest to see if we parametrize z{t) = x{t) + with 

x(to = 0) = Co, y(0) = 0, X, yreal. (3.3) 


Note that the symmetry of the contour under reflection in the real axis forces x to be an even 
function of t (so that x'(0) = 0), and y to be an odd function. We can rescale t to make y'{0) = 1. 
The expansion in eqn. (3.2) is then 

(x'(0)2 - y\0f + 2ix'{0)y'{0)) + ■ ■ ■ 

(3.4) 


F{zit)) ~ F(co) + 


^^(co)-q^.^T... 

Since all the derivatives of F{x) are real, the equation lmF{z{t)) = 0 is then satisfied to this 
order, and as F"(co) is positive, the function decreases with t. However, for conventional choices 
of /(x,(5o)) the contour does not continue parallel to the imaginary axis; to see where it does go, 
we need to consider the expansion to one higher order. 


F{z{t)) ~ F{co) - ^ 1 


r + • • • 


To lmF{z{t)) = 0 then requires 


x"(0) = 


F(3)(co) 


(3.5) 


(3.6) 


(3.7) 


3F"(co) 

Thus in the neighborhood of Cq, the contour has the form 

. X . F(^Hco) 2 

■ 

{F^^\co) is typically negative for the class of functions in which we are interested, as is necessary 
for this to be a useful truncation.) What is not so obvious — but turns out to be true — is that 
for our purposes, the contour C' described by this equation is essentially indistinguishable from 
the true contour of steepest descent. That is, in the region where the integrands of interest give 
the bulk of the contributions to the inverse Mellin transform, the two contours are extremely close, 
and so we incur almost no efficiency penalty in choosing the simplified contour given by eqn. (3.7). 
With this choice, our inverse Mellin transform can now be written 


1 

vr 


dt Re l^l - F{z{t)) . (3.8) 

As we shall see in later sections, this quadratic contour is a significant improvement over the 
linear contour chosen in ref. [7], not to mention the textbook contour. 
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4. Evaluating the Inverse Transform 


The next point we must consider is the evaluation of the inverse Mellin transform using our 
new contour (3.7). We want to choose a few points along the contour at which to evaluate the 
function in order to approximate the integral by a finite sum, 



(4.1) 


In order to find a ‘good’ set of points, we should in principle first find a set of functions in 
which to expand A possible (though not necessarily “optimal”) choice is a set of orthogonal 

polynomials with an astutely-chosen weight function. In fact, for integrating such a choice 

is optimal if we pick the weight function to be F{z) itself! We might therefore be tempted to 
construct a set of orthogonal polynomials with respect to this weight function. 

First, though, let us examine the behavior of the integrand along the contour near cq. Assume 
that /^{Qq) does not vary exponentially, so that it cannot eliminate the exponential fall-off expected 
from the x~^ factor. We then expect a behavior of the form F{co)e~^^*\ where the power series 
expansion of g{t) can be determined by matching coefficients with the power series expansion in the 
integrand of eqn. (3.8). (Note that the imaginary part of z'{t) in the prefactor does not contribute 
until order .) We then obtain 


9{t) 




(4.2) 


2 F(co) 

(As discussed in the previous section, both F"{cq) and F{cq) will be positive.) This suggests that 
we perform a change of variables u = whereupon the integral of eqn. (3.8) becomes 


1 / F(co) ^ 

7r\ 2F"{co) Jo y/u 


Re 


1 — i 


. E(3)(co) 2F{co) 
3F"(co) V F^co) 


y/u F Z 



2F(co) 

F"{co) 



= ^ I -^e “ Re [e“ (l - m 2 C 3 \/u) F [z {c 2 y/u))] 


du 


(4.3) 


27r 


where C 2 = y^2F(co)/F"(co) and C 3 = F^^^(co)/(3F"(co)). For small u, we expect the factor inside 
the brackets on the second line to vary slowly. What is again not so obvious, but turns out to be 
true for initial parton distributions of interest, is that the factor in brackets varies very little, and 
smoothly at that, in the entire region where the integral receives noticeable contributions. It is 
thus an excellent candidate for approximation by polynomials. The same statements continue to 
hold for the inverse Mellin transform of in which the factor inside the brackets in eqn. (4.3) 

is now multiplied by the evolution operator E^{as{Q'^),as{Qo)). 

The reader may also recognize the prefactor in front of the brackets as the weight function for 
the generalized Laguerre polynomials Ln ^^^^(u), and it is a generalized Gauss-Laguerre quadrature 



formula employing these polynomials which we should use for evaluating the integral in eqn. (4.3). 
This formula approximates an integral 


pOO 1 ^ 

/ h{u) ~ '^Wjh{uj ), (4.4) 

Jo Vu 

where the Uj are the zeros of Ln and the weights are given by standard formulm [11], 

(4.5) 


r(n + l/2) 
n! (n + 1)^ 


^n+l y^j) 


1 2 • 


5. Singlet Evolution 


The same considerations discussed in the previous two sections apply to singlet evolution, with 
some important differences. In the non-singlet case, the integrand in the inverse Mellin transform 
of f^iQo) is modified by a simple (and modest) multiplicative factor to obtain the integrand for 
f{x,Q‘^). In contrast, in the singlet case, the evolution operator has non-trivial matrix structure, 
and is not close to the identity matrix even for near Q^. This happens because the (S,^) basis 
is not an eigenbasis even for the lowest-order Altarelli-Parisi function. Accordingly, contours chosen 
according to either or will not be particularly good ones. What we want is a basis in which 
the evolution operator does not twist one direction into another as we move around in the complex 
plane. Such a basis is given by the eigenvectors of the lowest-order Altarelli-Parisi matrix, that is 
by Thus in the singlet sector, we could rewrite our Mellin integral as 


ic 


where 


Re —idzx ao)s^((5o) = 


'c 


Re —idzx ^{as,ao)Pls^{Q'^) + Re —idzx ^{as,aQ)P^s^{Q'^) , 


(5.1) 


'c 


s^(Q") = 


(5.2) 


'S-(Q2) 

g-{Q^) 

In each of these integrals, each of the components is now modified multiplicatively (up to 
0{as) corrections), like the non-singlet Mellin integrals I analyzed in previous sections. We might 
expect them to be treatable in the same fashion — contours chosen according to the components of 
P^s^{Qq). Indeed, following this approach we would choose four different contours for the different 
integrations above. There are, however, two complications we would confront. 

Analyticity assures us that it is legitimate to choose different contours for the PI and P^ 
integrals, but there is a subtlety: the integrands have branch cuts in the complex plane, and the 
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contours chosen according to section 3 may cross these branch cuts. So long as we use the same 
contour for both integrals, this is completely harmless, because all that happens upon crossing a 
branch cut is that interchange roles. This merely interchanges the two integrands, and the sum 
is unaffected. If we want to choose different contours, however, we should either shift the branch 
cuts so that neither integral crosses them (which is possible only some of the time), or else we 
must compute the integral around the branch cut. While the latter is possible, it will not allow an 
efficient evaluation of the integral. 

The other complication is that the components of P^s^{Q‘q) are no longer necessarily positive 
functions. As a result, they no longer necessarily have a minimum to the right of the rightmost 
pole. (Typically and will have nearly the same rightmost pole, as it controls the x —> 0 
behavior; the gluon distribution will be slightly steeper [1].) Indeed, in practice one combination 
— Pf s*((5q) — is roughly a sum of the two components, and hence is positive, and can be handled 
according to the prescriptions of the previous section; while the other combination is 

not, and in fact has no minimum to the right of the rightmost pole for certain choices of /^(Qg). 
(The contour of steepest descent heads straight into the pole.) 

Thus instead it is better to pick a slightly different approach to choosing the contour-determining 
function in this case. Instead of the parton distributions at the original scale Qq, take the distribu¬ 
tions evolved (using leading-order evolution) to another fixed scale as the contour-determining 
function F{z). The latter is of course just given by the evolution operator EIq{Q‘1, Qq) multiplied 
by the starting distributions s^(Qg). While in principle we might choose different contours for the 
two different components in the singlet sector, in practice it is more efficient overall to choose a sin¬ 
gle contour using the gluon component in the role of F{z). That is, for the purposes of determining 
the contour of integration, and the integration points along it, take 



(asiQDV^-'^^ , o. (Ps{Ql) 

V ao ) ^ n ao 



(5.3) 


for F{z) in the singlet sector. In this expression, Qi would be a scale intermediate between Qg and 
the highest scale to which we wish to evolve the parton distributions. 


6. Numerical Performance 

At how many points do we need to evaluate the integrands constructed according to the 
prescriptions in previous sections, in order to obtain an evolved parton distribution function to 
a given accuracy? While the answer will depend on the precise form of the initial distribution 
functions, and on the desired accuracy, we can obtain a very good idea by studying the numerics of 
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the evolution of the toy parton distribution set considered by Blumlein et al. [12], Those authors 
consider the MS massless four-flavor evolution from a reference scale Qq = 2 GeV, with = 
250 MeV. The initial distributions are taken to have the following form, 


Ql) = - xf, d^{x, Ql) = AdX~^-^{l - x)'^, 

S{x,Ql) = Asx~^-^{l-xy, c{x,Ql) = 0, (6.1) 

g{x,Ql) = Agx~^-^{1- xf . 


The sea S{x,Q‘^) = [T, — — dy] {x,Q‘^) includes the charm content, and is taken to carry 15% 

of the nucleon momentum at the input scale. This fixes As; the remainining Ai are fixed by the 
flavor and momentum sum rules. 

Take the desired accuracy to be the same 2 parts in 10^ sought by the above authors. For the 
non-singlet evolution (for example, the valence up or down densities) to Q = 10 GeV requires around 
30 points at small x, rising to around 80 points for x = 0.7, using Gauss-Legendre quadrature on 
the contour of ref. [7]. (Evolution to Q = 10 TeV requires about 10% more points at larger x.) 
In constrast, using contour derived in the present work requires only three or four points for all 
10“® < x < 0.9. In the singlet sector, the new contour typically requires four points, while the 
contour of ref. [7] requires roughly the same number of points as the non-singlet sector at small 
X, and somewhat more at larger x. Obtaining the toy parton set to the stated accuracy typically 
requires more points, because the charm distribution is a small difference of two larger numbers. 
The quadratic contour derived in this paper requires a lone additional point for most x values, 
and a total of seven points at the largest x value. (For the contour of ref. [7], the toy parton set 
typically requires 30 points at small x to over 120 at the largest x value.) 

The answers obtained using this new approach should be compared with those of the ‘truncated 
analytic solution,’ that is the lower half, of the table of ref. [12]. The following table illustrates 
the convergence of sample values using the new contour, with one, two, or four Gauss-Laguerre 
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points 


x 

/ 

xf{x) 

n = 1 

2 

4 

exact 

10“^ 

Uy 

frac. error 

0.0096683 

2.7-10-2 

0.0094433 

3.4-10-3 

0.0094111 
-2.7- 10-3 

0.0094114 

10-3 

Uy 

frac. error 

0.089636 
1.8 • 10-2 

0.088177 
1.0 - 10-3 

0.088077 
-1.0- 10-4 

0.088086 

0.1 

Uy 

frac. error 

0.47445 
3.8 • 10-3 

0.47291 
5.2 - 10-^ 

0.47267 
7.4- 10-^ 

0.47267 

0.3 

Uy 

frac. error 

0.31061 
8.6 • 10-3 

0.30793 
-1.2 - 10-^ 

0.30797 
4.1 - 10-3 

0.30797 

10-3 

9 

frac. error 

102.92 
4.9 • 10-2 

97.879 

-2.0-10-3 

98.080 
4.4- 10-^ 

98.080 

10-3 

9 

frac. error 

22.052 
4.5 • 10-2 

21.106 
-3.2 - lO-'^ 

21.112 
6.8- 10-3 

21.112 

0.1 

9 

frac. error 

1.4168 
1.2 • 10-3 

1.4152 
1.5 - 10-3 

1.4151 
1.5 - 10-^ 

1.4151 

0.3 

9 

frac. error 

0.18849 
5.0 • 10-3 

0.18754 

-3.8-10-3 

0.18755 
-2.2 - 10-^ 

0.18755 


7. Crossing Functions 


Crossing functions arise naturally in a crossing-symmetric formalism [13] for evaluating general 
next-to-leading order distributions in a collider environment. They give the change in the differential 
cross section as a colored final-state particle is crossed from the final state into the initial state. 
The crossing functions Ca^p{x, Q^) are factorization-scheme dependent. They can be expressed in 
terms of scheme-independent functions Aa^p and scheme-dependent functions Ba^p as follows, 


Ca^p{x,Q‘^) = 


N 


Aa^p{x,Q‘^)\n 




+ Ba^p{x,Q‘^) 


(7.1) 


where Smin is the boundary between the real contributions integrated analytically and those inte¬ 
grated numerically. A physical quantity will be independent of Smin, in the limit that Smin ^ 0. 
The A and B functions are convolutions of universal kernels with the parton distribution functions, 


Aa^p{x, = K^^f^ix) <S) fb^p{x, , 

Ba^p{x, = A:iLb(x) (g) fb^p{x, , 


(7.2) 
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(again with implicit summation over b). Explicit expressions for the kernels may be extracted 

from ref. [13], 


- .) + 2 (+ (1^ + .(1 - .) 


6iV 


(1 -x)+ 


X 


^q^q'i^) ^qq' ( ^ ^2 






-A 


K^^,ix) = - 1- 


J_^^ 1 + (1-Z)2 
iV2 


K; 


_)_ 2 ^^ _j_ 1)2(1 _ a;) ^ 



II 

C>1 

-Q 

<if(x) 

1 

“ m 




1 

Jp 



--jli(i-.)-5(i-.) + -(i+.=)(l^ 


1 

Jp 


1 + (1 - xf 


ln(l — x) + X 


(7.3) 


(Note that etc.) In this equation the ( )+ prescription is defined by 


{F{z))+ = hm^ ^(9(1 - z - P)F{z) - 6{l - z - P) F{y) dy^ 


(7.4) 


such that, 


= /‘&«khLhl)+5(i)i„g(i_^), 


(1 - z)+ 


1 - z 


c \ 1 - z /+ Jx ^-z 2 


(7.5) 


provided that g(z) is a function well behaved at z = 1. 
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The Mellin moments of the kernels are 


K. 




+ 


z{z - 1) {z + l){z + 2)y ’ 


= 6gg' ( 1 - 


q^q 




iV2 


+ 




2N \z + 2 z{z + l){z + 2)) ' 


kA,z ^ - I _ 




1 


+ 


z-1 z{z-l){z + l) 


<4"’' = + Ilf + «-(z + 1) - ?*(1))“ + V'(l) - <^'(2 + 1) 


3 

+ 


18 m 

2 

+ 


z{z-l) ' {z + l){z + 2)) V^(^ + l)) (^ + 1)2 + (^ + 2)2 


^q^q' ~ 


^B,MS _ _}_ 

q^g 2N 


1 - 


1 


TT^ 7 1 

- T + 


N^J [6 4 ' 2^2 2z(z +1)2 + 

+ 2z(zVl) (^(^) " ~ 


1 


+ 


z + 2 z{z + 1 )( 2 ; + 2) 
2 

+ 


(V^(l) -i;{z + 3)) 
1 

+ 


K. 


B,MS,z _ 


g^q 


= ^ l- 


Ar2 


+ 


z{z-\-2) z(z-\-l)(z-{-2) j 
2 


z-1 z{z-l){z + l) 


(V'(i) - V'(z + 1)) + 


+ 


{z + 1)2 z{z — 1) 


These moments allow us to write 


-)2\ _ t^B,z 




(7.6) 


(7.7) 


and then using eqn. (2.9), 




(7.8) 


We can evaluate the crossing function, given by the inverse Mellin transform of this moment, 
using the same approach described in previous sections, but with /^(Qo) replaced by 
{X = A, B). 

The conjugation identities mentioned above, along with the fact that K^^-, = 0, imply that 
the non-singlet and singlet sectors do not mix in eqn. (7.8), and that the kernel in the non-singlet 
sector is simply , while the 2x2 matrix in the singlet sector i 


IS 


/^X,2 J^X,Z 

I q^q q^g 
\ g^q g^g 


(7.9) 
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8. Comparison with x-Space Codes 


As mentioned in the introduction, the two main methods now commonly used to evolve parton 
distribution functions are the Mellin transform method [7,14] pursued in previous sections, or direct 
integration [15] of the differential equations for f{x,Q‘^). While formally equivalent to NLO, these 
techniques differ in the higher-order terms implicitly retained. As noted by Bliimlein et al. [5], 
these differences can have a substantial effect on the parton distributions, as much as several 
percent even at moderate x. One lesson from this comparison is that a determination of q;s(M^) to 
1% will require NNLO calculations; in the meantime, however, it is useful to be able to reproduce 
answers obtained through the x-space method as well. The purpose of this section is to show 
that solutions identical to those obtained by direct integration can be obtained using the Mellin 
transform technique, taking advantage of the same enhancements considered in previous sections. 

Let us first consider as itself; in terms of the QCD scale parameter A, it can be written [6] at 
NLO as 


asiQ'^) = 


1 




( 8 . 1 ) 


/3olnQ2/A2 /32 InQ^A^ 
let us replace it with an approximation (to NLO) for in terms of Ug at a different scale, instead 
of one in terms of Aqcd, 


as{Q‘^,Ql,ao) = 


ao(l +aQl3ohiQ‘^/Ql) 


[1 + ao(3ohiQ‘^/Qlf + ao/3r- (1 + ao/^r- + ao/3o InQ^/Qo) 


In 


1 + 


ao/3olnQ^/Qg 

A+aoPr) 


( 8 . 2 ) 

where j 3 r = /?i//9o) and the number of flavors is taken to be constant throughout the interval [Qo, Q]- 
This particular form of Ug emerges by performing one Newton-Raphson improvement to the 1-loop 
solution. We could recover the solution of the implicit (beta function) equation for Og, 


1 

agiQ'^) 


1 

asiQl) 


+ /loln [Q^Ql) 

Po 


asiQ‘^){(3o + PiagjQl) 
aa{Ql){/3o +Piag{Q‘^) 


(8.3) 


if we iterate the evolution'^, 

a'"' (g^ Ql ao) = (g^ Q?, (Q?, QI, ao)), 

4°'(g^go>«o) = ag{Q‘^,Ql,ao) . 


(8.4) 


where Qi lies in the interval [Qq, Q]; since the evolution is logarithmic, we would presumably choose 
Qi = \fQQo- Formally, we want to take the limit n ^ oo, so that 


<(g^gol«o) 


lim af\Q^,Ql,aQ). 

n —>-oo 


(8.5) 


t Iteration will result in convergence to the numerical solution of eqn. (2.6) only if the form of the approximation 
is suitable; for a slowly varying function whose true zero is near the initial approximation, a Newton-Raphson 
approximation is indeed suitable. 
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In practice, however, this is a completely pointless exercise, because the form (8.2) is already within 
one part in 10 “^ of the ‘exact’ answer for Q > Qq, so long as Qo > 1.5 GeV; in the other direction, 
for Qo = Mz, the fractional error is less than 2 • 10“^ for Q > 4 GeV. (The ‘exact’ answer refers to 
the numerical solution of the implicit NLO equation (8.3).) As Q grows for fixed Qo, the fractional 
error reaches a maximum for some Q'^, and then falls off; for Qo = 2 GeV, it is as noted less than 
1 part in 10“^, around Qm = 27 GeV. 

With a formula for as that reproduces a direct integration of the (3 function in hand, we can turn 
to the evolution equations for the parton distributions. Equation (2.7), for the Mellin transform of 
the evolved parton distribution, is identical to the original x-space equation. In the usual Mellin 
space approach, one performs the further truncations discussed in section 2, effectively dropping 
terms which are 0(af), to obtain the usual solution (2.11,2.13). I should stress that from a physical 
point of view, to next-to-leading order the truncated equations are no less valid than the original, 
and the solutions qualify no less as NLO solutions. What is of interest is the discrepancy between 
the two solutions: one or both has higher-order corrections at least as large as half the difference 
between them. In order to study this difference, we must solve the untruncated equation (2.7). 
Introduce an evolution operator just as in equation (2.9). In the non-singlet case, we can integrate 
the equation directly, obtaining 


In the case at hand. 


£'"‘(as,ao) = exp 

P^{a) 


P^{a) 


' ao 


/3(a) 


da = — 


LJ ao /^(®) 

Pq + aPf 


' ao 


a{Po + Pia] 


da 


da 


( 8 . 6 ) 


pz 

-in 


Po 


= -l^ln Li _ LT In 


ao 


-f’ 1 1 f Po + Pias \ 


(8.7) 


Pi \/3o + Piao) 


so that 


E^as,ao) = - 

\ao 


( 8 . 8 ) 


as\ f Po + Pias \ 

V/3o + Pio-o ) 

It is worth noting that this evolution operator is not much more expensive to evaluate on a computer 
than the conventional one, eqn. (2.11). 

Now, although we can solve the equation in closed form in the non-singlet case, we will not be 


able to do so in the singlet case, so it is useful to understand what alternative means we have of 
computing the operator in eqn. (8.8). Since it is an evolution operator, we can write it in the form 


E^{as,ao) = n^=iE^(aj,aj_i), 


with an = as, where (for example). 


(8.9) 


aj — ao P (as ao)j/n . 


( 8 . 10 ) 
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Equation (8.9) is, of course, exact; but can we take advantage of it to approximate E^{aj, aj-i)7 
We would then hope to recover the full evolution operator by taking the n ^ oo limit, 

E^{as,ao) = lim aj_i). (8-11) 

In fact, we can, but there is a subtlety concerning the expansion parameter. Define 5a = {as—ao)/n, 
the total evolution length for each step in eqn. (8.9), and expand E^(aj,aj-i) in 6a, 

E^ (nj,(ij_i) = l-t-ei(nj)(In-t-e2(nj)((^n)^“t“*** (8.12) 


We then find 


E^{as,aQ) = 1 + 6a'^ei{aj) + {5af ^ ei(ajJei(ajJ H-h (5a)”ei(ai) • • • 61 ( 0 ^) 

j=l = i 

Ji#i2 

n n 

+ 62(0^) + ( 5 a)^ 62(«ji)e2(aj2) + ''' + (<J«)^"'e2(ai) • • • 62(0^) + • • • 

j = l Jl02 = l 

il#32 

/ n \- 

= 1 + (5a ei(aj) 1 + terms with fewer sums than powers of 6a . 

\ / 

(8.13) 

If the Ci are well-behaved functions (as they are in our case), roughly speaking each complete sum 
over an index produces a factor of n compensating the 1/n implicit in 6a, so that when we take the 
limit the terms with fewer sums than powers will vanish because they are suppressed by powers of 
n. Thus in the limit we obtain 


E (a^, uq) 


lim 

n —>-oo 


1 + 6a J^ei(a, 
i=i 


exp 


/■““ dlnE^(a,ao) 

L — 


(8.14) 


which using (2.7) is of course equivalent to eqn. (8.6). So we learn that expansion in 6a is legiti¬ 
mate, and that the higher-order terms don’t matter in the limit (they would of course accelerate 
convergence if they were present). This is neither surprising nor subtle; the subtlety comes when 
we consider in addition expanding in aj, which is also a small parameter. 


Gi{a) — Ci^o + e^qa -|- -!-••• 


(8.15) 


and truncating at order m, to obtain ei{a). Running through the above argument, we find that 
no matter how large the number of segments n gets, there is always an error of 0(a™'"'“^) in the 
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estimate of E^{as,ao) 


E {CLs ; 0^0 ) 


lim 1 + (5a > ei( 

r ).—I * ^ 


exp 


da 


LJ ao 


j = l 

d\nE^{a, ao) 
da 


+ 0(a 


m+l\ 


(8.16) 


Unfortunately, the standard evolution operator (2.11) does involve precisely these sorts of trunca¬ 
tions, so we can’t use it to recover the missing higher-order terms. We can obtain a kernel that 
will work by exponentiating only the leading-order term of the integral (8.7), and leaving the rest 
expanded, 

/a 

E^{as,ao) = (^fJ 


Pi V/^o + Picio ) 


(8.17) 


This is equivalent to (formally) expanding in which is in general not a small parameter; the 
expansion only works because the accompanying logarithm is small. We may thus expect that the 
speed of convergence varies as we move around the complex plane, and the question then arises of 
how large n has to be in equation (8.9) before we converge to the ‘exact’ answer given by eqn. (8.8). 
The answer, somewhat amusingly, is that for the non-singlet parton distribution functions one 
encounters in practice, a real-world value of a^, and evolution over the range from Qq = 2 GeV to 
Q = 15 GeV, n = 1 suffices (except at the smallest values of x, for which we need n = 3). That is, 
using eqn. (8.17) already gives the limiting answer, to within 2 parts in 10“^. More generally, n = 2 
limits the error to 6 parts in 10'^ for Q < 10 TeV. 

In the singlet sector, equation (2.7) has the solution 


E (a^, ao) 


T^exp 



P^ja) 

p{a) 


da 


(8.18) 


where the a-ordering operator T^, the analog of the usual path-ordering operator, orders matrices 
P^{a) according to their distance from ao, putting matrices of arguments further away from ao 
further to the left. It appears because in general P^ matrices at different points a, a' do not 
commute. 

Equation (8.18) is still a formal expression, and one must make further approximations to 
obtain an explicit expression. We can make use of eqn. (8.9), in which for singlet evolution we must 
also order the matrices, 

E^{as,ao) = E^{an,an-i)---E^{ai,ao). (8.19) 


Let us take a 5a sufficiently small that we can formally expand NLO terms in the a-ordered 
exponential (8.18) to first order, without assuming that the same is true of LO terms; we can then 
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^ ^ f r r* da 

E (as,ao) = r„ exp / — 

I L A^O Jao ® 

^ r \ n da 


- ^ _ (ji_\ 
/3(<.) UftoJ 


( 8 . 20 ) 


Evaluating the time-ordered product gives us our final expression for the singlet evolution operator, 
E (ug, ag ) 

/^ \ -A'-//3o f 1 


P_ - In 


/^o + /dlQ-s 

/3o + /?! ®0 


-P_P1P_^ 


/3o “ -|- A1 


(ai-a;)//3o 


Fi . 1,1- 


/3o ’ Po ’ Po ' 


( 8 . 21 ) 


/ Ai - A^ Ai - A^ /3i \ 

+ (+ ^ “) 

In this equation, 2 E 1 is the hypergeometric function; while Pi is substantially larger than f3o, the 
argument of the hypergeometric function will still be much smaller (in absolute magnitude) than 
one, so that we can evaluate it using its series expansion. 


, u N r(c r(n-ha r n-f 6 ^ 

2 F 1 (a, b] c; z) = > -————- 2 : 

r a)r(6 ^ E n + c n! 


( 8 . 22 ) 


To reproduce the toy parton as discussed in section 6, we would need a rather large number of 
subintervals n in eqn. (8.19) for evolution at small and large x. It is possible to reduce this number 
by resumming the terms in Pi that are proportional to Eg. To do so, define 


Then 


M^_=Tr P_Pi , /i; = TT P+P 


Pi = /ilP_ + P-P 1 P+ + + P+PiP- , 


(8.23) 


(8.24) 


and 

E^(as,ag) = 

PPo + PiaA-^-^^^ 


P_ -P_PiP+ 


1 / /3i \ 

13,-K + Ai (' + a“V 


/3o + / 1^ /3o — -h \ /3o / 

^ ^ Ai-Al A^-Ai /3i 

- 2 T 1 i + M^-Ml, 1-^-3—;2-^-3— 

V« 0 / V Po Po Po 

( z z K-^P /3l \ 1 1 

-ao 2U (1 + 1 - 2 - —-^a„) 


+ (+ ^ -) 


(8.25) 
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With this form, 20 subintervals suffice at x ~ 10“^, 10 or fewer at intermediate x, and only at the 
largest x ~ 0.7 do we still require a large number (> 100) of subintervals. 

Using equations (8.2,8.8,8.25,8.19) yields the evolution as would be given by a direct integration 
of the untruncated evolution equation (2.1) along with a direct integration of the beta function (2.6). 
One might have thought that this would be the evolution as computed by various x-space programs; 
but it isn’t. The reason is that these programs do not appear to integrate the beta function 
numerically, but rather use the approximate solution (8.1). This provides yet a third inequivalent 
version of NLO evolution. It isn’t difficult to reproduce it in the Mellin approach, however; avoid 
changing variables to a^, and instead integrate eqn. (2.5) with respect to Q^. In the non-singlet 
sector, we obtain 


E^{Q\Ql) = 

-PoV/3o 

L 


exp 




1 


pz 


In Ln In L 


_£]_pz 

2/34 1 


1 -|- 2 In To 


To L 
1 -|- 2 In L 




+ 


A. 

27(31 


■PI 


L2 

2 -|- 6 In Lg -|- 9 In^ Lg 


(8.26) 


2 -|- 6 In L -|- 9 In^ L 
Z3 


where the hat on E indicates that its arguments are momentum scales instead of running couplings, 
and where Li = InQf/A^. We can obtain the analog to equation (8.17) by formally expanding the 
exponential in the Pf, but leaving the power prefactor unexpanded. In the singlet sector, the 
first-order expansion in eqn. (8.20) yields 


E^{Q\Ql) = r 


P_ +P_ 


-A^_//3o 


1 


02^1 


1 

^ \Lo 

PI 


- 

PA 


In Ln In L 


Lt] 


L 


Pi pz 

A ' 


1 -|- 2 In Lg 1 -|- 2 In L 


+ P- 


+ 


1 


27/3, 


-Pi 


0 


2 -|- 6 In Lg -|- 9 In^ Ln 2 -|- 6 In L -|- 9 In^ L 




L3 


-P( 


- 2 


+ 2 


l_/3o(/3o +~-^+) ^ V-^o E 

Pi pz( 1 

/32(2/3o + Ai-A-+)2 1 

Pi 


rL 


- 2 


Pi 


P( 


El 


P_ 
In Lg 


L2 


InL 


El 


1 

- 


PAPo + Ai - A^) 


pz 

2^1 


L2'^ 
In Lg In L 


Poi'^Po + Al — A^) ^ 

I 

Pl{3(3o + Xt-Xl) ' 

PI 


El 


L2 ^ 

In^ Lg In^ L 




L3 




L3 


+2 


P( 


/3g3(3/3g + Ai-A^)3 1 V^o E^ 




+ 


+ (+--) 


(8.27) 

where = L/Lg and = (A^ — X(_)/(3q. This form of the evolution operator should again be 
used in combination with eqn. (8.19); it does require a rather large number of subintervals n if we 
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want an answer accurate to 2 parts in 10"^. It is thus better to add in the terms arising from an 
expansion of the time-ordered exponential to second order, 


5E^^{Q\Ql) 


C’’l,2,3 = ± 






\ 2 p2 

’ ^(73 ? ^ (T2 


)Pc 


0-3 


(8.28) 


The (somewhat lengthy) formula for T 2 is given in appendix II; it should be kept in mind that most 
of the computation need only be done once, not anew for each value of . 

With the addition of the terms in r 2 , most values of x require n = 3 or 4 in eqn. (8.19), though 
at x = 0.7, 30 subintervals are required to provide an accuracy of 2 • 10“^ for the toy parton set of 
ref. [12]. 

Using these evolution operators, along with eqn. (8.1), I find that I indeed reproduce the 
results of the ‘direct solution,’ that is the upper half of table 1 in ref. [12], mostly within their 
quoted errors. Curiously enough, the values of the evolved parton densities are quite close to what 
would be obtained from the ‘exact’ evolution described earlier in this section, typically within a 
few parts in 10^. Note, however, that the running coupling values are somewhat different; with 
Qo = 2 GeV, the two values for the running coupling differ by about 1% at Q = 10 GeV. 


9. Conclusions 

In this paper, I have derived a quadratic contour for calculating the evolution of parton distri¬ 
bution functions within the Mellin transform method, and demonstrated its superiority over other 
techniques in the literature. I have also shown how to reproduce the results obtained within the 
‘x-space’ method using a modified evolution operator. In addition to the application discussed 
here, the method may also be used within a general framework for extracting parton distribution 
functions from collider data [16]. 

I thank M. Peskin and A. Vogt for helpful comments. 


Appendix I. Analytic Continuation of a Dilogarithm Integral 

In order to evaluate the evolution operators along our chosen contours in the complex plane, 
we must be able to evaluate the moments of the leading- and next-to-leading order anomalous 
dimensions at essentially arbitrary points in the complex plane. As discussed in ref. [7], all functions 
appearing in the Mellin moments of these anomalous dimensions, save one, have expressions in terms 
of elementary functions or derivatives of the P function. The latter can be calculated everywhere 
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in the complex plane via well-known techniques. The one exception is S{z), for which GRV give 
the following expression, 


Siz) = --a 3 ) + v 


{'ipiz + 1 ) - C( 2 ) 


(tpiiz + l)/2) -'i(j{z/2)) + / dxx 




where r/ = ±1, ■0 is the usual logarithmic derivative of the gamma function, and Li 2 is the diloga¬ 
rithm. As it is not completely clear that the expansion given by the authors of ref. [7] for the last 
term is valid for our purposes, I give an alternate one here. 


Define 


Mi(z + 1) = [ 


Using ^ = (1 + x)x^ ^ — x^ we find the following recurrence relation for Mi, 

Mi{z + 1) =—Mi{z) + f dx x^~‘^Li2ix) 

Jo 

= Mi{z — 1 ) — / dxx^~^{l — x)Li 2 {x). 

Jo 

The other integrals in this equation may be performed by integration by parts. Let 

M2{z) = [ dx x^~^Li2{x) = -I- ^ — [ dxx^~‘^ln{l — x) 

Jo z -1 z -1 Jo 


z-1 z -1 da ^^0 Jo 


dxx"" ^(1 - x)“ 




M2{z) — M2{z — 1) = — f dxx^ ^(1 —x)Li2(x) 

Jo 


(z - l){z - 2) V (z - 2)2 (z - 1)2 J (z- 1)3 • 

(1.5) 

On the other hand, we can also write down an ‘asymptotic’ expansion for Mi; to do so, 
repeatedly rewrite 

1 1 — X 1 

1 + x ~ 2{l+x) 2 ’ ^ 


to obtain 


Mi(z) = - f dxx^ ^Li 2 (x)- f dxx^ ^{1 — x) 

2 Jo 2 Jg 


Li2(x) 

1 + x 


-^2 ^ / dxx^ ^(1 — x)'^Li2(x) -|- 2 dxx^ ^(l — x)^ 


Li2(x) 


Af i / ■ \ j 

V2-^’V(-l)' (-J ) M2 (z + 0+ 2-^-1 / dxx^-^{l-x)^+^- 

7=0 1=0 -^0 ^ 


Li2(x) 
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Since in the z ^ oo limit, the integral is dominated by the region x ~ 1, the last term goes 
as and we can drop it if we are interested in the asymptotic expansion only through 

order + 1. In fact, because of the factor of 2~^ in front, dropping this term is a reasonable 
approximation even for modest z. We can re-expand the resulting expression in terms of l/z, but 
this yields a rather poor representation. Instead we can rewrite the remainder term 


)-7V-l 


Jo 1 + X 

Li 2 (l) / dxx^-\l-x)^+^ + 2 
Jo 


■,-N-2 


-N-1 


dxx^ '^{1 — x) 


fLi2{x) xLi2(l) 


= 2 


-N-2 


C(2) 


r(z)r(A^ + 2) 
r(z + A^ + 2) 




dxx^ ^(1 — x) 


1 + X 2 

N+i /^ Li2(x) _ xLi2(l) ^ 

1 + X 2 


( 1 . 8 ) 

The second integral goes as \nz/z^'^^. For Rez sufficiently large, we simply drop the second term; 
for other z, we can use the recurrence relation (1.3) [backwards] to shift z into this range. It’s 
worth using the recurrence relation (1.3) a few times explicitly, since we will certainly be interested 
in values near z ~ 1 or 2 for which the approximation with reasonable N won’t be enough; we can 
then simplify the transcendental functions to minimize the number of such function evaluations we 
need to perform. 


Appendix II. Second-Order Terms for the x-Space Singlet Evolution Operator 


To present the (somewhat lengthy) formula for T 2 , used in equation (8.28), define the following 
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functions 


d .. =c+— _ 

Po 

^11 = ’ 

Po 

Ao = , 

Po 

Pi 

Poi = ^PoPPi , 

Po 

^00 = , 

Po 

= Al - dl;23Ao 
pP = Pn - (ii;i 2 Poi 
nP = ^01 - di,23Poo 
Pw = PlO - dl;12P00 
s- =\-1 Po- 


As in section 8, tl = L/Lq. 



Then 


(Ai, A 2 , A 3 , p) = - ^h 23 lnLor/ p Pl^ 0 ^dl ;12 - Aod 2;13 - ^d 2;13 + 2 Pood 2 ; 13 ) 

PqLq \ / 


+ {pit^dr,,2 - Pit^d2;i3 - 4“^d?;i2 + Piod%33 + Pl>i^d: 


2;13 ~ 2PoOC?2;13 


Pi II QI’l / p(a) 1 1 I p(^) J p(“)^ J p J p J j2 

- /34 r 3 - I -^11 “1;23«2;12 + «1;12«2;23 ~ Ml «1;23“3;13 ~ Ml«2;23«3;13 ~ M0«l;23«2;12 

PoM ^ 

— -Poi(^l;12rf2;23 + -foi<^3;13'i2;23 + 2-PlO'^l;23(^3;13 + 2-Poi<^2;23'^3;13^ 

- p r3 (-f’ll^^l;23<^2;12 + -Pll^'^l;12'^2;23 “ -Pll'^3;13'i2;23 “ -^11 ^I^l;23'^3;13 “ -Pllrf2;23'i3;13 

Po-^o V 

+ -Poi'i2;23'^3;13 + 2T’lO'il;23d3;13 + 2Poi'^2;23C^3;13) 

- aeri ^ {‘^Plld2-,23d\.i2 + 2-Plld2;12'i2;23 “ 2Pii(i4;i3d2.23 + -^11 ^I^l;23<^3;12 “ -Pl0<^l;23(^3;12 

PO-^0 ^ 

+ P[^^d 

1;12<^3;23 -flll^4;13<i3;23 -foi<^l;12<^3;23 +'^01<^4;13'^3;23 -^11 ^'^l;23l^4;13 

— 4Piid2;23<^4;13 ~ -f’ll'^3;23(^4;13 + 2Tbl'^3;23'^4;13 + 3-PlO'il;23'^4;13 + 3-PoiI^3;23'^4;13 

Oy02 ^1 

- 4)6 r 4 (^-^I1^2;12'^2;23 + -f’ll^'^l;23(^3;12 + -Pi?'^l;12(^3;23 “ -Plll^4;13'^3;23 “ 2Piid2.23d4;13 

Po-^0 ^ 

— -fll(j^3;23*^4;13 4" '^0ll^3;23^4;13 “ ^11 ^I^l;23'^4;13 “ 4:i^il(i2;23<^4;13 

— -PllC?3;23<^4;13 + 2 ^ 01 ^ 3 . 23 (^ 4.13 + 3Pio(il;23(^4;13 + 3Poi(^3;23(^4;13 


^8p5 ^ ('^2;23^3;12 + '^2;23'i3;12 + '^2;12'^3;23 + I^2;12'^3;23 '^5;13'i3;23 '^2;23'^5;13 

-PiiC?3;23(^6;13/3iI’ 


2i^ 3;23'^5;13 3(i2;23'3^5;13 3(i3;23'^5;13) + 
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4PiiAirA /p p 


llA^l'L /'^2 43 I 42 ^3 43 j2 j2 ^3 r,^2 j3 04 44 04 44 

/08r5 M2;23“3;12 “T “2;12“'3;23 “3;23“5;13 “'2;23“5;13 ^“3;23l''5;13 'J“2;23n5;l3 'J0'3;23«'5;13 

Po -^0 

4Pii(i3;23 IiiLqAiI’l M 43 I 42 42 42 42 04 43 f;44 \ 

4 'pToJe M3;23«3;12 + “3;12“3;23 ^3;23^6-,13 >Jn3;23'46;13 D®6;13; 

4Pii(i3;23AlI’L /43 42 42 43 0 4 44 fjjS \ 

“I ^ 10^6 M3;12‘l3;23 “3;23‘l6;13 'J“3;23«6;13 '4“6;13; 

(il;12(il;23l’L / „(a) „(a) j \ , ^’or^(^l;12l^l;23 ln , i^lo^(il;12dl;23 ln 

Mio + LI3IL„ 

^ 111 ^0'’t _j_ 3i.23<l2:l2^I^ 2/?l fil0l(l;23<(2;12 ^ 111 


+ 


LAo'Lo ' I2p4^^ 

2AlPii^(il;23C^i;12l’L 2AlPlorfl;23 ln ^orfi;12l’L 2P*^“^di;23 ln ^rfi;12/3l I’l 








2T’lodl;23 In L In Lodli2Pl'^L 2P^^“'’(il;23rf3-12/3lI’L 2Pio(il;23 In Lo(i3.i2AlI’L 

+- -r Q /nfi T T- O /nfi t- -+ 




L^fdlLo 


L^PlLo 
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2 / 3 i-PiV^ 1 ; 12 '^ 2;23 ‘^PiPmdi-i2d2]2‘i^^L\nL qv^ 2l3iP[i di-i2d\2ir^ 

' T /~iA T O T /^A T O ' 


mhi 


mil 


mil 


2(5iPoidi-^i2 lnZ/(i|.23r ^ 4Pii(i2;i2<^2;23 InLInL q/?! ^Piid2-,2‘i Pod\i20i'^i 


Lf3^L^ 


0-^0 


L^^lLl 


0-^0 
<5? 


L^PlLl 


4Piiii2;i2lnL<i|23/5i^L 4Piid|32(^2;23/5i 4Pii(j2 ;23 In L In Lcdl.^Plr^L 


L^(5tLl 


L^fdtLl 


L^(5lLl 


4Pii hlLdl.22,dl.i2[dlr^L 4Pii(i2;23 ln-Lorf3;12/5l’’L 4Pii(i|23'^3;12/3l’’L 

“I” T O /n« T O —^ +- T- O /nC T- O —^ ' 


L^(5lLl 


L^fdlLl 


2 Pii^(il ;12 lnLorf3;23/5l’’L 2 Poidl ;12 In L In Lo(i3.23/3ir^^ 2P3^3^<il;i2li3;23/3l^L 


L 3 / 38 L 2 

yib). 




Lfd^Lf. 




2-Poi'^ 1;12 lnLd|23/3ir^^ 4Pii(i2;i2 lnLlnLo(i§.23/3ir^^ 4Pii lnLod|i2'^3;23/^i^L 

“I” T- /n« T- Q -^ + - -r n /nQ -r “3 -+ - 


L^qLq 


4 Pii(i 2 ;i 2 lnLd|. 23 / 3 ir^" ^ 4 Pii(i|. 32 d 3 ; 23 /?i^L ^Pn InLInLo(ii.i 2 '^i; 23 / 5 i’’L 

TTT^sTo I” 


L 2 / 38 L 3 

:{ 3 j 2 


L^PlLl 

34 


L 2 / 38 L 3 


L 2 / 38 L 3 


L 3 / 310 L 3 


4Pii lnLoc^3;i2'^3;23/5i’’L ^Pn InLd§.i2d3;23/5i ’’l 2Pii<i3;23 In L(ii.i2/3ir^^In^^Q 


P 3 / 310 L 3 

^2 

di-,23d2-,i3 In Pr/ 


L^Pl^Ll 


L^Pl^Ll 


Pio + Pqi^ ~ 2 Pood 2 ;l 3 ) + 


'^ 1 ; 23 <^ 2 ; 13 ^L / p(a) 


L 2/32 


PlT - Pi 00 ! 2;13 - P^fd2-13 + 2 Pood 2 - 


^ 2;13 


— -^2 4-^ ^PiV^1;23 + Plld2-23 ~ 2Pio(il;23C^3;13 ~ 2Poi(i2;23C?3;13 ~ Poi'i2;23) 


P 3 / 3 , 
“^Pids-Asrp 


^ (^Pn^ di-23d3-13 + Plld2-23d3-,13 + Pll'i2;23 “ Poi<^3;13C?2;23 “ 2Piodi;23rf3;i3 “ 2Poid2;23<^3;13 


+ 


P^/^ePo 

P^^Po 

L’dlLl 

2Pll(i3;12(^2;23/^l^L 

2Pllli3;12 In Po(i3;23/34L 

2Plld3;12di;23/?^fln'i 

L3/3§L§ 

L3/3‘"Li> 

i3,3l«L3 

Pl^l^dl;12d3;23/?frf 14 Po ^ 

Poi(^l;12'^3;23 In P/3^r^^ In^Po 

2Pll(i2;12'^3;23 In L^fr pllpLo 


Lpm 


Lpm 


Plld3;12d3;23/3f?’L In^Pln^Po C^4;13/3 ??’l In^P ^ D ^ , D ^ ^ , 

- 56 - l-rioai ;23 + p 3 lM 3 ; 23 j + 


L 3 / 310^3 


L 4 / 3 g 


^ 0-^0 

2 Plld 3 ; 23 di;i 2 / 3 l?'L In^Po 


^Piid^-isPfrpln^L 2Pnd3-23de-i3Pfrpln^L ‘^Piid^i2d3-23(dif'L 

- ^5^8 - (^2;23 + d3;23) + - ^6/310 - (^3;23 + ^de-w) “ 

2 Pn/ 3 frj 4 n 3 Lo 


PlLl 

(dtLt 


P 6/330 

('^2;23<^3;12 + '^2;12<i3;23 ~ <i2;23'i5;13 ~ 'i3;23<^5;13) ~ 


L^Pl^Ll 

2 Pll(i 3 ; 234 l 2 / 3 ^Lln 2 Lo 


P 3 / 310 L 3 


(Pi 0 ' 3 ^ 1 ; 23 '^ 3;12 + Poi'^l; 12 (^ 3;23 “ P’l 0 C^l; 23 '^ 4;13 “ Poi' 3 ? 3 ; 23 ^^ 4 ; 13 ) 


26 



+ 


2d4;13 In ( ^(a) 

L4/36 


-Pil'’'^l;23C?4;13 + 4Pll(i2;23rf4;13 + -Plld3;23rf4;13 + 2^11^2.23 + -Pll'^3;23 “ 2Poirf4;13C?3;23 

— -Poi<^3;23 “ 3Piodi;23rf4;i3 “ 3PoiC^3;23'^4;13) 


H-^4|^6 ^ ^2Piid4;i3(i2;23 + -Pllrf4;13C?3;23 + -Pll'^3;23 “ -foi<^4;13C?3;23 + -PlV^l;23'^4;13 + 4Pll(i2;23C^4;13 


+ -Pll'^3;23C?4;13 “ 2Poi<^3;23'^4;13 “ 3Pio(il;23'^4;13 “ 3PoiC?3;23'^4;13) 


4Pii(i5;i3 lnL(]frl^ 


(c? 5;13<^2;23 + 2d5;13C?3;23 + ^3;23 + 3(i2;23'^5;13 + 3li3;23C^5;13) 


4-Plld5;13/34r^ ^ 2 j2 I j3 I QJ j2 , qj ^2 ^ 

-^^^^8- (“5;i3d2.23 + 2a5;i3a3.23 + a3;23 + 3a2;23«5;13 + 3a3;23«5;13) 


+ 


+ 


4Pii(i3;23 lnLd^.i3/3^r^=* 


L 6/310 


(3li3;23'^6;13 + 'i3;23 + 6'^6;13) + 


2/3id3;i3rinn'L 

L3/34 


(Piodi ;23 + -foi<^2 ;23) 


Pmdi-23d2-,i3'r L 4:Piid^.2‘id^.i^f3fr£ 


ijl 


L^(5f 


+ 


0 

H-- ^-Pll^<il;23 + 4Pll(i2;23 + Plld3;23 ~ 3Piolil;23'^4;13 ~ 3Poi(^3;23'^4;13 ~ 2Poi^^3;23) 

2Piid5;i3/3?rfln2L 


L 6/310 


(3(i3;23C?6;13 + '3?3;23 + 6d6;13) 


L5/3i 


0 


+ 


+ 


+ 


‘^Piid3-23de-i3(dir £ \t? L 

2(3ir^j^\nLo 

/^o-^o 

(dtLt 


(3(^2; 23 <^ 5 ; 13 + 3(i3;23C?5;13 + <3^123 + ‘^dl. 23 ) 

(3d3;23'i6;13 + <^3523 + 6 d 6 ;i 3 ) H- ^’^2 r ^2 -^ ('^ 1;12 “ <^2;13) 

Po-^0 


{Plodl;23d2;12 + -Poi'il;12'^2;23 ~ -Pl0^^1;23<^3;13 ~ -foi<^2;23^^3;13) 

^4Pii(i2;12'i2;23 + -^11 ^(^1;23'^3;12 + -Pll^rfl;12rf3;23 “ -Pll^rfl;23C?4;13 “ 4Pii(i2;23C?4;13 


— Plld3-23di;13 — 2Plodl;23'^3;12 “ 2Poi(^l;12'^3;23 + 2-Poi(^4;13'^3;23 + 3-PlO'^l;23C?4;13 + 3-Poi'^3;23(^4;i3^ 
^5 -- ('i3;23'i2;12 + ^3;12C?2;23 “ '^5;13C^2;23 + 2li2;23'^3;12 + 2(i2;12'i3;23 “ 2(i5;i3(i3.23 

Po-^0 ^ 

oj j2 oj j2 \ , -Plld3;23/3iri In'^Lo ^ ^ 

- 3d2;23«5;13 “ 3d3;23«5;13 j 4- ^iQ^g - (®3;12 “ a6;13) 


+ 


2Piici3;23/3frpln"Lo 


(2d3;23'i3;l2 + '2^3;12 + d3-,12d\.23 — 'i6;13'^3;23 “ 3(i3;23C?6;13 “ 6 ^ 6 . 43 ) 


+ ^^^^^3;23^^1^L -^0 (^ 3 , 42 ^ 3.23 - (i3;23C?6;13 + ^ 3:12 “ 2d6;13) ' 

Po -^0 


(II.2) 
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